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ABSTRACT 

Two dimensional modeling of collisionless shocks has been of tremendous importance in understanding 
the physics of the non-linear evolution, momentum transfer and particle acceleration, but current computer 
capacities have now reached a point where three dimensional modeling is becoming feasible. We present the 
first three dimensional model of a fully developed and relaxed relativistic ion-electron shock, and analyze and 
compare it to similar 2D models. Quantitative and qualitative differences are found with respect to the two- 
dimensional models. The shock jump conditions are of course different, because of the extra degree of freedom, 
but in addition it is found that strong parallel electric fields develop at the shock interface, the level of magnetic 
field energy is lower, and the non-thermal particle distribution is shallower with a powerlaw index of '^2.2. 

Subject headings: acceleration of particles — instabilities — magnetic fields — plasmas — shock waves 
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In collisionless shocks the mean free path of individual 
particles is much larger than than the characteristic scales 
in the shock structure, and effective collisions and pressure 
support are not mediated in particle-particle interactions, but 
by collective forces. Collisionless shocks are ubiquitous in 
the universe, and happen on a range of scales from the bow- 
shock of the solar wind at the earth, over solar corona, su- 
pernova remnants to highly relativistic shocks at gamma-ray 
bursts and active galactic nuclei. Given that only collective 
forces act as an effective collision mechanism, instabilities 
can easily arise in the plasma. It has been shown over the 
last ten years how a range of Weibel-li ke filamentation insta- 
bilities op erate at t he shock interface ("Kazimura et al."1998'; 
Medvedev & Loeb 1999; Bret et al. 2004; Frederiksen et al. 



2004 



ISilva et al. 2003 ; Nishikawa et al. 2003; Ke shet et al.| 



2009), which can lead to magnetic field generation and par- 
ticle acceleration. A breakthrough for the study and under- 
standing of collisionless shocks has been the application of 
ab initio particle-in-cell simulations, where the formation and 
evolution beyond the linear phase of the instabilities of the 
shocks can be studied with almost no assumptions. Until now, 
the largest scale computer experiments of both magnetized 
dShimada & Hoshino| [20041 |Saito & Sakai| [20041 'Hededal' 
& Nishikawa||2005[ |Spitkovsky||20Q5t T 
2009D and unmagnetized (Haugb0lle 7 



Sironi & S pitkovsky 
2005; 'Hededal 2005; 
^piSbvsky 2005; Silva 2006; Kato 2007; Chang et al. 2008; 
Spitkovsky 2008a b; Martins et al. 2009) collisionless shocks 
nave mostly been performed in two dimensions, while three 
dimensional simulations of unmagnetized shocks have been 
too costly to scale to large enough sizes. Some of the first 
studies of electron-positron and electron-ion shocks were per- 
formed in 3D, but only the very early linear and quasi-linear 
stages of the shock formation could be followed (Frederiksen, 
|erSl[2QQ4l ISilva et al.|[2Q03l [Nishikawa et"aLl|2003| ). Three' 
dimensional studies of pair plasmas have been done showing 
the full development of the shock ramp, and recovering the 
correct jump conditions fHaugb0lle" 2OO5[ [Spitkovsky "2005 '; 
[Nishikawa et al. 2009), but only limited electron-ion simula- 
tions have been made (Spitkovsky 2008a). 

The development of long term large scale two dimensional 
simulations have helped tremendously in the understanding 

[haugboel @ nbi. dk| 



of collisionless shocks, both in quantitative and qualitative 
terms, but given current computational resources, and the 
technical quality of the Particle-In-Cell codes, we now have 
the possibility of using fully three dimensional modeling to 
properly account for the full dynamics. In this paper we 
present the first three dimensional simulation of a relativis- 
tic unmagnetized collisionless electron-ion shock, where the 
simulation is followed long enough to establish the correct 
jump conditions, create a fully thermalized and developed 
downstream region, and follow the emergence of a power- 
law distributed population of high energy ions and electrons 
downstream of the shock. We compare the results to similar 
2D simulations, and show how the different dimensionality 
impacts on the formation and evolution of the shock. 

1. SIMULATION SETUP 

We have used the massively parallel 3D particle-in-cell 
Photon-Plasma code, developed at the Niels Bohr Institute 
(Haugb0lle 2005; Hededal 2005), to simulate collisionless 
electron-ion shocks. For the three dimensional ion-electron 
run we used a 2^^ order field solver and TSC interpolation of 
particles, while for all other runs we used a newly developed 
6* order field solver and cubic interpolation of p articles , giv- 
ing effectively a 50% higher resolution (see e.g. [PCCP| ). The 
simulation domain is periodic perpendicular to the streaming 
direction, while the lower inflow boundary is open for out- 
going particles and has absorbing field boundaries and the 
opposite, upper boundary is a perfectly reflecting wall. We 
use a combination of a mo ving injector, to launch the shock 
( [Sironi & Spitkovsky 2009|), and a moving window, to follow 
the shock evolution ( Tzeng et al.|1 996). Initially, only a small 
slice in front of the wall is populated with a streaming plasma, 
and a shock is launched when the reflected fields and particles 
collide with the inflowing plasma. New particles are added to 
the plasma in front of the shock with the speed of light un- 
til the full box is populated with a streaming plasma. At the 
lower boundary outgoing particles are allowed to escape the 
box, and electromagnetic fields are damped while the shock 
transition propagates through the box and the downstream re- 
gion grows. When the shock transition reaches the center 
of the box, a moving frame is applied which continuously 
moves the simulation domain to keep the shock transition at 
the center. Initially the plasma is unmagnetized, and hence 
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electron skin depths 

1/2 1/2 

1 . — Upper panel: From top to bottom the ion density, the magnetic field density e^ , electric field density e^ , and the electric field projected along the 
magnetic field lines E^^g at oupe t = 2250. The fields are normalized to the kinetic energy density of the upstream bulk flow, as for example eg = B^ /[2n{mi + 
Me) (F- 1)]. Lower panel: From top to bottom the ion density (normalized to 1/3 upstream), the magnetic field density, and the electric field density at ujpe t = 1500 
for a 2D coUisionless shock. 



nearly transparent for the particles. The particles reflected 
by the wall creates an artificial counter- streaming population 
of particles upstream of the shock transition region. While 
in 2D simulations the domain can be made large enough for 
this initial spike of particles to diminish in density and be- 
come insignificant, the required domain sizes are currently 
prohibitively expensive in 3D. But with the combination of 
the open lower boundary and the moving frame technique we 
can nevertheless follow the evolution of the shock until it has 
settled to a steady state, with a thermalized population down- 
stream of the shock, and an in- streaming population mixed 
with particles reflected from the shock upstream. All runs 
have been done with an upstream Lorentz factor of F = 15. 
The 3D simulation was done with 250x250x7000 cells, 6 par- 
ticles per species per cell, a mass ratio mi/rrie of 16, and con- 
tained up to 16 billion particles. To properly resolve the dy- 
namics both up- and downstream of the shock we resolve the 
electron skin depth Se = {Airnee^ /rUeTy^^ with 14 grid cells 
upstream and 7 cells downstream of the shock, corresponding 
to 56 and 23 cells per ion skin depth up- and downstream of 
the shock. For the time stepping we used a Courant condi- 
tion of 0.4, checking both light crossing time in a grid cell, 
and the local plasma frequency. To check the impact the box 
size had on the evolution we ran 2D simulations with initially 
250x7000 cells and then wider and longer boxes with up to 
1000x14000 cells. Additionally, to probe a larger dynamic 
range in 3D we also performed a pair-plasma simulation with 
375x375x3750 cells and 7.5 cells / skin depth, for a volume 



of 50x50x500 skin depths. 

2. SHOCK STRUCTURE AND EVOLUTION 

Figure [T| shows the structure of the evolved 3-D shock at 
ujpe t = 2250. Upstream of the shock the particles reflected 
at the shock interfaces interact with the incoming particles, 
generating a two-stream instability where the ions collect into 
current channels and are pinched by the self-generated mag- 
netic field. The electrons Debye- shield the ions while oscil- 
lating in the strong transverse electric field, giving effective 
electron heating and momentum transfer from ions to elec- 
trons. The current channels slowly merge while loosing mo- 
mentum, and break up at the shock interface. The behavior 
of the upstream medium in the 3D simulation is in qualita- 
tive agreement with the 2D simulations, though the current 
channels in 3D seem to merge less than in 2D. This may be 
due to the fact that current channels in 2D (being ID curves) 
necessarily cross, while in 3D that is not the case (compare 
e.g. densities in fig. [T]). The average eg and the level of elec- 
tron to ion momentum are similar in 2D and 3D (see figure 
t, while the average eg with a maximum at 8% is slightly 
wer in 3D. Eventually the highly intermittent electromag- 
netic fields in the upstream current channels have transferred 
sufficient bulk flow energy to heat and transverse momentum 
to make the channels diffuse and slow down. A transition to 
a downstream shocked medium occurs over approximately 20 
ion skin depths, equivalent to about one Larmor radius (given 
that the ions and electrons have approximately the same en- 
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Fig. 2. — Average energy distribution across the shock at cjpet = 2250. To 
avoid biasing only particles with positive velocities are used in calculating 
the average kinetic energy per particle, and the magnetic field density is nor- 
malized to the in- streaming bulk kinetic energy. 
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Fig. 3. — Particle distribution function sampled downstream of the shock 
above the dashed line marked on figure|5]at oupe t = 2250. To the left are elec- 
trons and to the right ions. The red dasned line is a relativistic Maxwellian, 
the purple line is a powerlaw with a low and high-7 cut-off, and the blue line 
is the combined model. 

ergy, their Larmor radii are similar). While the width of the 
shock transition is nearly the same in 2D and 3D, unique for 
3D is the occurrence of large electric fields parallel with the 
magnetic field (see upper panel in figure [T]) in the transition 
region, coupled with eg = 1 locally, giving an effective mech- 
anism for ion and electron acceleration across the shock tran- 
sition. A similar phenomena has been observed with satellites 
in the collisionless shocks in the aurora and the Earth's mag- 
netosphere, where parallel electric fields together with strong 
particle acceleration are observed at the shock transition (see 
e. g. Ergun et al. 2001 , 2009 j . We also note that parallel elec- 
tric fields play a keyr ole in reconnection of magnetic fields 
( [Schindler et al.||1988| ), such as the reconfiguration that hap- 
pens here, going frorti the mostly transverse current driven 
axial magnetic field upstream of the shock to the turbulent 
flux ropes seen downstream of the shock. In the downstream 
region the plasma is very nearly neutral, with little density 
variation, but with a high level of magnetic turbulence. Our 
current simulation box for the ion-electron shock is not wide 
enough to allow for the largest magnetic structures, but using 
a 3D pair plasma simulation with a larger transverse volume 
we have observed how closed flux ropes are formed and are 
ad vected downwards from the shock, similar to what was seen 
by |Spitkovs^ ( [2QQ5] ). 



3. PARTICLE ACCELERATION AND DISTRIBUTION 

In 2D models of collisionless shocks it has been found that 
particles are sl owly acce lerated b y scattering off the filaments 
( [Hededal et al 2004; S pitkovsky 2008b ; Martins et al. 2009). 
A few "lucky" particles cross back and forth over the shock 
interfaces a number of times, and this Fermi-like accelera- 
tion process slowly builds up a power-law tail of high energy 
particles, due to the quasi constant probability that a single 



high-energy particle is r eflected in the stron g transverse elec- 
tric field near the shock (Martins et al. 2009| ). Even though the 
acceleration mechanism is not identical to the normal Fermi 
process, the resulting power-law index is in good agreement 
with theoretical predictions, and is approximately a =23 - 
2.6, where f(p) ex 7"" at high energies. In our 3D simulation 
we find the same emergence of a high energy tail distribution, 
on top of a relativistic Maxwellian downstream of the shock. 
We model it as 



/(V7) = Ai(v7)^exp (-m,-,7/r) +A27"'' 
[1 +exp(-(7^/„-7)/A^/„)]"^ 

[l-\-QXp(-(j-Jjnax)/^max)r^ 



(1) 



where A/ are normalizations, T is the temperature, a is the 
powerlaw slope, and 7/ & A/ are the locations and widths of 
the cut-offs. We impose the low energy cut-off so the power- 
law makes a smooth match to the Maxwellian, while the upper 
cut-off is time-dependent and grows with time. We find that 
in the 3D simulation the slope is shallo wer with o^ = 2.1 
2.3, m atching the theoretical prediction of |Keshet & Waxman 
( |2005| ), for a F = 15 relativistic shock of atheo = 2.1. 
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Fig. 4. — Density profile of the 2D and 3D runs. In the upper panel the 
density profile at cupe t = 1600 is shown for a 14000x1000 cell domain, and a 
7000x250 cell domain. In the lower panel the density profile for the 3D run a 
the times oopt = 1500 (red line) and oupe t = 2250 (blue line) is shown. 



4. CONVERGENCE 

The density and momentum distribution are some of the 
lowest order criteria to check when modeling a collisionless 
shock. To assess the impact of our relatively limited simula- 
tion domain we compared the different 2D runs, finding that 
the limited box size has only a minor impact on the jump 
conditions and on the evolution of the filaments. Analyzing 
the 2D runs we find that a box with 250x7000 cells contains 
a shock with velocity Vsh = 0.42c, and a downstream to up- 
stream density ratio of rid/riu = 3.24, while a box with twice 
the length, independent of the transverse size, more faithfully 
reproduces the analytic jump conditions. With Vgh = 0.47c and 
a density jump of 3.12 it is in percent precision agreement 
with the analytic expectation for a relativistic gas Ud/riu = 
7ad/[7ad - 1] + l/[r(7ad " 1)], where 7ad is the adiabatic in- 
dex of the gas, which is 4/3 (3/2) for a 3D (2D) relativistic 
gas. This is also seen in the 3D case, where we find Vsh = 0.27 
and Ud/riu = 4.62, where analytically one expects Vsh = 0.31 
and Hd/riu = 4.2 (see figure |4]). Because we launch the shock 
reflecting cold streaming particles on a wall it takes some time 
until a proper, thermalized downstream region is created. The 
electromagnetic fields at and near the shock interface have to 
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Fig. 5. — Evolution of the ions in phase space. The PDF is sampled to 
the right of the dotted line. Notice that it is only after uopt ~ 2000 that the 
downstream part of the shock is completely thermalized and at rest. At earlier 
times the impact of the wall is still significant, with the return current of 
reflected particles apparent in the camel shaped pdf. 

build up to a sufficient level to scatter the particles, and the 
shock interface has to be far enough away from the wall at 
the upper boundary, so that upstream particles are scattered 
by the fields and thermalize thoroughly before they possibly 
reach the wall. This convergence can be monitored by look- 
ing at the v^7 momentum near the wall. A camel shaped PDF 
signals that proper pressure support in the downstream region 
has still not been established. We find (see figure Bl that the 
simulation has to run to uOpe t = 2000 or Upt t = 500 before 
a proper equilibrium is established, even though the proper 
jump conditions are already established at Upet = 1000 (see 
figure |4]). Without the moving frame approach it would have 
required a box with a least 28000 cells in the streaming direc- 
tion to properly establish and thermalize the shock. 



5. CONCLUSIONS 

In this Letter we have studied for the first time the long time 
behavior of a fully developed 3D collisionless ion-electron 
shock. We find that the development of the shock structure 
is similar to that of an electron-ion shock in a 2D model, but 
there are both qualitative and quantitative differences of im- 
portance when making quantitative predictions from shock 
models for observations. Apart from the dimensionality, 
which gives rise to different shock jump conditions, the extra 
degree of freedom changes the shock in a number of ways: 1) 
The current channels upstream of the shock become more sta- 
ble. 2) Very strong parallel electric fields along the field lines 
are created at the shock interface giving a new and effective 
avenue for particle acceleration, as compared to 2D models. 
3) The level of magnetic field energy is lower near the shock 
interface 4) A power-law tail in the PDF downstream of the 
shock emerges at late times. In agreement with theory the 
power-law index is shallower in the 3D shock (~2.2) com- 
pared t o the 2D model (^2.5 ). 

In T rier Frederiksen et al.| (|2010) it was found that the ra- 
diation emitted from the Weibel instability in 2D and 3D is 
qualitatively different, and analogously we expect that the dif- 
ferences presented here in the physical development of a 3D 
shock compared to a 2D shock will give rise to significant 
differences in the emitted radiation. This is a topic for future 
studies. 
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